Project2

Gaussian random field with application of INLA

Part I Multivariate normal distribution

Let \(\boldsymbol{x} = (x_1, \dots, x_n), n = 100\) be multivariate normal distributed with \(E(x_i) = 0, Var(x_i) = 1\), and \(Corr(x_i, x_j) = e^{-0.1|i - j|}\)

1. Compute and image the covariance matrix \(\boldsymbol{\Sigma}\) of \(\boldsymbol{x}\)

2. Find the lower Cholesky factor \(\boldsymbol{L}\), such that \(\boldsymbol{L}\boldsymbol{L}^T = \boldsymbol{\Sigma}\), of this covariance matrix, and image.

3. Sample \(\boldsymbol{x} = \boldsymbol{Lz}\), where \(\boldsymbol{z}\) is a length n random vector of independent standard normal variables. Plot the sample.

4. Find the precision matrix \(\boldsymbol{Q}\) of the covariance matrix, and compute the lower Cholesky factor \(\boldsymbol{L}_Q\), such that \(\boldsymbol{L}_Q\boldsymbol{L}_Q^T = \boldsymbol{Q}\), of this matrix. Image these matrices and compare them to the images obtained in a) and b)

5. Sample \(\boldsymbol{x}\) by solving \(\boldsymbol{L}_Q^T\boldsymbol{x} = \boldsymbol{z}\), where \(\boldsymbol{z}\) is a length n random vector of independen standard normal variables. Plot the sample.

6. Permute the ordering of variables in \(\boldsymbol{x}\), and redo the exercises.

---

Solution to Part I

a)

Given that \(\boldsymbol{\Sigma} = e^{-0.1|i - j|}\). The covariance matrix can be expressed as follows: \[\begin{equation*} \Sigma = \begin{pmatrix} 1 & e^{- 0.1 h_{12}} &\dots & e^{- 0.1 h_{1n}} \\ e^{- 0.1 h_{21}} & 1 & \dots & e^{- 0.1 h_{2n}} \\ \vdots & \vdots & \ddots & \vdots \\ e^{- 0.1 h_{n1}} & e^{- 0.1 h_{n2}} & \dots & 1 \end{pmatrix} \end{equation*}\]

### b) According to the cholesky decomposition rule, \(\boldsymbol{L}\) is the lower triangular matrix for \(\boldsymbol{\Sigma}\), it can be easily computed from R using L = chol(Sigma). It is then plotted as below.

c)

Sample using \(\boldsymbol{x} = \boldsymbol{L}\boldsymbol{z}\) transforms the zero-mean, standard normal random variales to the random variables with the desired covariance matrix.

d)

The precision matrix \(\boldsymbol{Q}\) is the inverse of the covariance matrix \(\boldsymbol{\Sigma}\), it is computed using Q = solve(Sigma) in R. The three matrices are thereby depicted as follows. Since the covariance matrix is not singular, given that the it belongs to the Matern family, thus it is analytically guaranteed to have positive definite property. Therefore, both precision matrix and the lower triangular precision matrix exist.

e)

Similarly, the expected random samples can be generated using the inversion of the above formula, thus \(\boldsymbol{L}_Q^T\boldsymbol{x} = \boldsymbol{z}\)

f)

Permute \(\boldsymbol{x}\) to make randomise the ordering of the grid, the asscociated covariance matrix can be thereby modified in a sparse way.

Part II Gaussian random fields and Kriging

The purpose of this computer exercise is to give an introduction to parameter estima- tion and kriging for Gaussian random field models for spatial data.

We assume the following observation model on the unit square:

\[\begin{equation*} y(\boldsymbol{s}_j) = x(\boldsymbol{s}_j) + \epsilon_j, \ \ \ j = 1, \dots, N, \end{equation*}\]

where \(\epsilon_j ~ \mathop{N}(0, \tau^2)\) are independent measurement noise terms. Further, consider a Matérn covariance function for the Gaussian random field \(\boldsymbol{x(s)}\):

\[\begin{equation*} Cov(x(\boldsymbol{s}_i), x(\boldsymbol{s}_j)) = \Sigma_{i, j} = \sigma^2(1 + \phi h)\exp(-\phi h), \end{equation*}\]

where \(h\) denotes the Euclidean distance between the two sites \(\boldsymbol{s}_i\) and \(\boldsymbol{s}_j\).

We assume the mean increases with east and north coordinates as follows: \(\mu_j = \alpha((s_{j1} - 0.5) + (s_{j2} - 0.5))\), for site \(\boldsymbol{s}_j = (s_{j1}, s_{j2})\) on the unit square.

---

2.1 Simulation

Simulate \(N = 200\) random sites in the unit square and plot them. Form the covariance matrix using \(\sigma = 1, \phi = 10, \tau = 0.05\). Take its Cholesky decomposition and simulate dependent zero-mean Gaussian data variables, then add the mean using \(\alpha = 1\). Plot your observations.

---

2.2 Paramter estimation

We will now use the simulated data to estimate the model parameters \(\alpha, \sigma^2, \tau^2, \phi\) using maximum likelihood estimation. Iterate between the update for the mean parameter, and updating the covariance parameters. Monitor the likelihood function at each step of the algorithm to check convergence.

## [1] "0.895765013504915  , iter no is  0"
## [1] "0.304910507162466  , iter no is  1"
## [1] "0.270287466981051  , iter no is  2"
## [1] "0.222351840261664  , iter no is  3"
## [1] "0.202516701628354  , iter no is  4"

## 
## Estimated sigma is  1.125 ; True sigma is  1 
## Estimated phi is  10.54 ; True phi is  10 
## Estimated tau is  0.01016 ; True tau is  0.05 
## Estimated beta1 is  -0.93 ; True beta1 is  -1 
## Estimated beta2 is  1.34 ; True beta2 is  1 
## Estimated beta3 is  0.92 ; True beta3 is  1NULL

---

2.3 Kriging

We will now use the estimated model parameters to perform kriging prediction. Predict variables \(x(s)\), where predictions sites lie on a regular grid of size 25x25 for the unit square. Visualize the Kriging surface and the prediction standard error. Compare with the true field.

## [1] "The prediction error is  1.93020697440398"

---

Part III